# 1 "butyM.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "butyM.f90"
!/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!*  Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.    *
!*  Copyleft (!c) 2005-2011 Elseif Laboratory.                      *
!*  elseifkk_users.sf.net                                           *
!*                                                                  *
!*  All Rights Reversed.                                            *
!*                                                                  *
!*                                                                  *
!*  LICENSE AGREEMENT:                                              *
!*                                                                  *
!*  You agree that we don't agree to anything.                      *
!*                                                                  *
!*                                                                  *
!*  DISCLAIMER OF WARRANTIES:                                       *
!*                                                                  *
!*  This software is provided to you "AS IS," and the authors and   *
!*  contributors disclaim any warranty or liability obligations to  *
!*  you of any kind, implied or express, including but not limited  *
!*  to the warranties of merchantability, fitness for a particular  *
!*  purpose, and non-infringement.                                  *
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
!
! bitmap image utility; Mother version.
!
! TTD:
! 1) catch allocation failure
! 2) error code and error handling routines
!
! note:
! 1) some routines use one dimenstional accessing
! to 3D array, which may cause some problems
! when compiling with boundary check option
! turned on.
!
# 45 "butyM.f90"
module butyR8
# 70 "butyM.f90"
!
  use butyD
  implicit none
  integer*4,private::DEBUG=0
contains
  !------------------------------------------------------!
  subroutine median_filter_r8(n,img,m)
    ! TTD: 1) treat image boundaries correctly.
    ! note: 1) real version of sort cannot handle negatives.
    use kuty ! for sorting
    use muty
    integer*4,intent(in)::n(3)
    real(8),intent(inout)::img(:,:,:)
    integer*4,intent(in)::m(3)
    integer*4 ix,iy,iz
    integer*4 ilo(3),iup(3)
    real(8),allocatable::b(:)
    real*8,allocatable::out(:,:,:)
    integer*4 k,mx,ktop
    integer*4 iix1,iix2,iiy1,iiy2,iiz1,iiz2
    integer*4 iix,iiy,iiz

    real*8 pmin,pmax
    ! because sort for real cannot process negatives correctly
    ! call zerocut(product(n),img(1,1,1))



    call get_minmax_img_r8(product(n),img(1,1,1),pmin,pmax)

    if(pmin<0) then
       img=img-pmin
    end if

    ktop=(2*m(1)+1)*(2*m(2)+1)*(2*m(3)+1)
    allocate(b(ktop))
    allocate(out(n(1),n(2),n(3))) ! <<< vulnerable
    mx=m(1)*2+1
    ilo=1+m
    iup=n-m
    do iz=ilo(3),iup(3)
       iiz1=iz-m(3)
       iiz2=iz+m(3)
       do iy=ilo(2),iup(2)
          iiy1=iy-m(2)
          iiy2=iy+m(2)
          do ix=ilo(1),iup(1)
             iix1=ix-m(1)
             k=1
             do iiz=iiz1,iiz2
                do iiy=iiy1,iiy2
                   call cparray(b(k),img(iix1,iiy,iiz),mx)
                   k=k+mx
                end do
             end do
             call sort(ktop,b)
             out(ix,iy,iz)=b(ktop/2)
          end do
       end do
    end do
    call mask(1,n(1), 1,n(2), 1,m(3))
    call mask(1,n(1), 1,n(2), n(3)-m(3)+1,n(3))
    call mask(1,n(1), 1,m(2), m(3)+1,n(3)-m(3))
    call mask(1,n(1), n(2)-m(2)+1,n(2), m(3)+1,n(3)-m(3))
    call mask(1,m(1), m(2)+1,n(2)-m(2), m(3)+1,n(3)-m(3))
    call mask(n(1)-m(1)+1,n(1), m(2)+1,n(2)-m(2), m(3)+1,n(3)-m(3))
    img=out
    deallocate(out)
    deallocate(b)

    if(pmin<0) then
       img=img+pmin
    endif

  contains
    subroutine mask(xlo,xup,ylo,yup,zlo,zup)
      integer*4,intent(in)::xlo,xup,ylo,yup,zlo,zup
      do iz=zlo,zup
         iiz1=max(iz-m(3),1)
         iiz2=min(iz+m(3),n(3))
         do iy=ylo,yup
            iiy1=max(iy-m(2),1)
            iiy2=min(iy+m(2),n(2))
            do ix=xlo,xup
               iix1=max(ix-m(1),1)
               iix2=min(ix+m(1),n(1))
               k=1
               do iiz=iiz1,iiz2
                  do iiy=iiy1,iiy2
                     do iix=iix1,iix2
                        b(k)=img(iix,iiy,iiz)
                        k=k+1
                     end do
                  end do
               end do
               if(k<3) cycle
               call sort(k,b)
               out(ix,iy,iz)=b(k/2)
            end do
         end do
      end do
    end subroutine mask
  end subroutine
  !------------------------------------------------------!
  subroutine mask_img_r8(lo,up,img,lo_val,up_val)
    real*8,intent(in)::lo,up
    real(8),intent(inout)::img(:,:,:)
    real*8,intent(in)::lo_val,up_val
    where(img<=lo)
       img=lo_val
    elsewhere(img>=up)
       img=up_val
    end where
  end subroutine
  !------------------------------------------------------!
  subroutine get_avevar_img_r8(n,img,plo,pup,ave,var)
    integer*4,intent(in)::n(3)
    real(8),intent(in)::img(:,:,:)
    real*8,intent(in)::plo,pup
    real*8,intent(out)::ave,var
    integer*4 i,m
    real*8 s,s2,z
    s=0
    s2=0
    m=0
    do i=1,n(1)*n(2)*n(3)
       z=img(i,1,1)
       if(z<plo.or.z>pup) cycle
       s=s+z
       s2=s2+z**2
       m=m+1
    end do
    if(m/=0) then
       ave=s/dble(m)
       var=dsqrt(s2/dble(m)-ave**2)
    else
       ave=HUGE(img(1,1,1))
       var=ave
    endif
  end subroutine
  !------------------------------------------------------!
  subroutine get_2dhist_r8(n,img1,img2,plo,pup,nh,h12)
    integer*4,intent(in)::n(3)
    real(8),intent(in)::img1(:,:,:)
    real(8),intent(in)::img2(:,:,:)
    real*8,intent(in)::plo,pup
    integer*4,intent(in)::nh
    real*8,intent(out)::h12(nh,nh)
    integer*4 i,ik1,ik2
    real*8 a,b,k1,k2,fk1,fk2
    real*8 x1,x2
    a=(dble(nh-1))/(pup-plo)
    b=(plo*dble(nh)-pup)/(plo-pup)
    h12=0
    do i=1,product(n)
       x1=img1(i,1,1)
       if(x1<plo.or.x1>pup) cycle
       x2=img2(i,1,1)
       if(x2<plo.or.x2>pup) cycle
       k1=a*x1+b
       ik1=int(k1)
       if(ik1<=0.or.ik1>nh) cycle
       k2=a*x2+b
       ik2=int(k2)
       if(ik2<=0.or.ik2>nh) cycle
       fk1=k1-dble(ik1)
       fk2=k2-dble(ik2)
       h12(ik1,ik2)=h12(ik1,ik2)+(1-fk1)*(1-fk2)
       if(fk1/=0.and.ik1/=nh) h12(ik1+1,ik2)=h12(ik1+1,ik2)+fk1*(1-fk2)
       if(fk2==0) cycle
       if(ik2/=nh) h12(ik1,ik2+1)=h12(ik1,ik2+1)+(1-fk1)*fk2
       if(fk1==0) cycle
       if(ik1/=nh.and.ik2/=nh) h12(ik1+1,ik2+1)=h12(ik1+1,ik2+1)+fk1*fk2
    end do
  end subroutine
  !------------------------------------------------------!
  subroutine get_hist_r8(n,img,plo,pup,nh,h)
    integer*4,intent(in)::n(3)
    real(8),intent(in)::img(:,:,:)
    real*8,intent(in)::plo,pup
    integer*4,intent(in)::nh
    real*8,intent(out)::h(nh)
    integer*4 i,ik
    real*8 a,b,k,fk
    real*8 x
    a=(dble(nh-1))/(pup-plo)
    b=(plo*dble(nh)-pup)/(plo-pup)
    h=0
    do i=1,product(n)
       x=img(i,1,1)
       if(x<=plo.or.x>=pup) cycle
       k=a*x+b
       ik=int(k)
       fk=k-dble(ik)
       if(ik<=0.or.ik>nh) cycle
       h(ik)=h(ik)+(1-fk)
       if(fk/=0.and.ik==nh) cycle
       h(ik+1)=h(ik+1)+fk
    end do
  end subroutine
  !------------------------------------------------------!
  subroutine remaphist_img_r8(img,plo_in,pup_in,plo_out,pup_out)
    real(8),intent(inout)::img(:,:,:)
    real*8,intent(in)::plo_in,pup_in
    real*8,intent(in)::plo_out,pup_out
    real*8 a,b
    a=(pup_out-plo_out)/(pup_in-plo_in)
    b=(plo_in*pup_out-pup_in*plo_out)/(plo_in-pup_in)
    img=img*a+b
  end subroutine
  !------------------------------------------------------!
  subroutine eqhist_img_r8(n,img,plo,pup,lo_val,up_val,nh)
    ! histogram equalization
    integer*4,intent(in)::n(3)
    real(8),intent(inout)::img(:,:,:)
    real*8,intent(in)::plo,pup
    real*8,intent(in)::lo_val,up_val
    integer*4,intent(in)::nh
    real*8,allocatable::h(:)
    integer*4 ik,i
    real*8 a,b,fk,p,aa,bb
    allocate(h(nh))
    call get_hist_r8(n,img,plo,pup,nh,h)
    do i=2,nh
       h(i)=h(i)+h(i-1)
    end do
    a=(dble(nh-1))/(pup-plo)
    b=(plo*dble(nh)-pup)/(plo-pup)
    aa=(pup-plo)/(h(nh)-h(1))
    bb=(h(1)*pup-h(nh)*plo)/(h(1)-h(nh))
    do i=1,product(n)
       p=img(i,1,1)*a+b
       ik=int(p)
       fk=p-dble(ik)
       if(ik>0.and.ik<=nh)then
          if(fk/=0.and.ik==nh) cycle
          img(i,1,1)=(h(ik)*(1-fk)+h(ik+1)*fk)*aa+bb
       elseif(ik<=0) then
          img(i,1,1)=lo_val
       else
          img(i,1,1)=up_val
       endif
    end do
    deallocate(h)
  end subroutine
  !------------------------------------------------------!
  subroutine rebin_img_r8(in,out)
    ! The origin of the image or the center of the image
    ! remains the same after rebinning.
    type(image_data),intent(in)::in
    type(image_data),intent(out)::out
    integer*4 ix,iy,iz
    integer*4 ixx,iyy,izz
    real*8 s
    real*8 up(3),lo(3)
    integer*4 iup(3),ilo(3)
    real*8 fup(3),flo(3)
    real*8 iin_off(3),out_off(3)
    real*8 norm
    out%n=idnint(dble(in%n)/out%d*in%d)
    call alloc_image_data(out,out%n,(64))
    if(iand(out%stat,IMG_ALLOCED)==0) return
    iin_off=dble(in%n)/2
    out_off=out%d*dble(out%n)/2
    norm=dble(product(out%n))/dble(product(in%n))
    do iz=1,out%n(3)
       lo(3)=((iz-1)*out%d(3)-out_off(3))/in%d(3)+iin_off(3)
       up(3)=(iz*out%d(3)-out_off(3))/in%d(3)+iin_off(3)
       if(lo(3)<0) lo(3)=0
       if(up(3)>in%n(3)) up(3)=in%n(3)
       ilo(3)=int(lo(3))
       flo(3)=lo(3)-dble(ilo(3))
       if(flo(3)/=0) then
          ilo(3)=ilo(3)+1
          flo(3)=1-flo(3)
       endif
       iup(3)=int(up(3))
       fup(3)=up(3)-dble(iup(3))
       do iy=1,out%n(2)
          lo(2)=((iy-1)*out%d(2)-out_off(2))/in%d(2)+iin_off(2)
          up(2)=(iy*out%d(2)-out_off(2))/in%d(2)+iin_off(2)
          if(lo(2)<0) lo(2)=0
          if(up(2)>in%n(2)) up(2)=in%n(2)
          ilo(2)=int(lo(2))
          flo(2)=lo(2)-dble(ilo(2))
          if(flo(2)/=0) then
             ilo(2)=ilo(2)+1
             flo(2)=1-flo(2)
          endif
          iup(2)=int(up(2))
          fup(2)=up(2)-dble(iup(2))
          do ix=1,out%n(1)
             lo(1)=((ix-1)*out%d(1)-out_off(1))/in%d(1)+iin_off(1)
             up(1)=(ix*out%d(1)-out_off(1))/in%d(1)+iin_off(1)
             if(lo(1)<0) lo(1)=0
             if(up(1)>in%n(1)) up(1)=in%n(1)
             ilo(1)=int(lo(1))
             flo(1)=lo(1)-dble(ilo(1))
             if(flo(1)/=0) then
                ilo(1)=ilo(1)+1
                flo(1)=1-flo(1)
             endif
             iup(1)=int(up(1))
             fup(1)=up(1)-dble(iup(1))
             s=0
             ! for bottom plane
             if(ilo(3)/=0)then
                if(flo(2)/=0) then
                   if(flo(1)/=0) &
                        s=s+flo(3)*flo(2)*flo(1)*in%img_r8(ilo(1),ilo(2),ilo(3))
                   do ixx=ilo(1)+1,iup(1)
                      s=s+flo(3)*flo(2)*in%img_r8(ixx,ilo(2),ilo(3))
                   end do
                   if(fup(1)/=0) &
                        s=s+flo(3)*flo(2)*fup(1)&
                        *in%img_r8(iup(1)+1,ilo(2),ilo(3))
                endif
                !
                do iyy=ilo(2)+1,iup(2)
                   if(flo(1)/=0) &
                        s=s+flo(3)*flo(1)*in%img_r8(ilo(1),iyy,ilo(3))
                   do ixx=ilo(1)+1,iup(1)
                      s=s+flo(3)*in%img_r8(ixx,iyy,ilo(3))
                   end do
                   if(fup(1)/=0) &
                        s=s+flo(3)*fup(1)*in%img_r8(iup(1)+1,iyy,ilo(3))
                end do
                !
                if(fup(2)/=0) then
                   if(flo(1)/=0) &
                        s=s+flo(3)*fup(2)*flo(1)&
                        *in%img_r8(ilo(1),iup(2)+1,ilo(3))
                   do ixx=ilo(1)+1,iup(1)
                      s=s+flo(3)*fup(2)*in%img_r8(ixx,iup(2)+1,ilo(3))
                   end do
                   if(fup(1)/=0) &
                        s=s+flo(3)*fup(2)*fup(1)&
                        *in%img_r8(iup(1)+1,iup(2)+1,ilo(3))
                end if
             endif
             ! for middle planes
             do izz=ilo(3)+1,iup(3)
                !
                if(flo(2)/=0) then
                   if(flo(1)/=0) &
                        s=s+flo(2)*flo(1)*in%img_r8(ilo(1),ilo(2),izz)
                   do ixx=ilo(1)+1,iup(1)
                      s=s+flo(2)*in%img_r8(ixx,ilo(2),izz)
                   end do
                   if(fup(1)/=0) &
                        s=s+flo(2)*fup(1)*in%img_r8(iup(1)+1,ilo(2),izz)
                endif
                !
                do iyy=ilo(2)+1,iup(2)
                   if(flo(1)/=0) &
                        s=s+flo(1)*in%img_r8(ilo(1),iyy,izz)
                   do ixx=ilo(1)+1,iup(1)
                      s=s+in%img_r8(ixx,iyy,izz)
                   end do
                   if(fup(1)/=0) &
                        s=s+fup(1)*in%img_r8(iup(1)+1,iyy,izz)
                end do
                !
                if(fup(2)/=0) then
                   if(flo(1)/=0) &
                        s=s+fup(2)*flo(1)*in%img_r8(ilo(1),iup(2)+1,izz)
                   do ixx=ilo(1)+1,iup(1)
                      s=s+fup(2)*in%img_r8(ixx,iup(2)+1,izz)
                   end do
                   if(fup(1)/=0) &
                        s=s+fup(2)*fup(1)*in%img_r8(iup(1)+1,iup(2)+1,izz)
                end if
                !
             end do
             ! for top plane
             if(fup(3)/=0) then
                if(flo(2)/=0) then
                   if(flo(1)/=0) &
                        s=s+fup(3)*flo(2)*flo(1)*in%img_r8(ilo(1),ilo(2),iup(3)+1)
                   do ixx=ilo(1)+1,iup(1)
                      s=s+fup(3)*flo(2)*in%img_r8(ixx,ilo(2),iup(3)+1)
                   end do
                   if(fup(1)/=0) &
                        s=s+fup(3)*flo(2)*fup(1)&
                        *in%img_r8(iup(1)+1,ilo(2),iup(3)+1)
                endif
                !
                do iyy=ilo(2)+1,iup(2)
                   if(flo(1)/=0) &
                        s=s+fup(3)*flo(1)*in%img_r8(ilo(1),iyy,iup(3)+1)
                   do ixx=ilo(1)+1,iup(1)
                      s=s+fup(3)*in%img_r8(ixx,iyy,iup(3)+1)
                   end do
                   if(fup(1)/=0) &
                        s=s+fup(3)*fup(1)*in%img_r8(iup(1)+1,iyy,iup(3)+1)
                end do
                !
                if(fup(2)/=0) then
                   if(flo(1)/=0) &
                        s=s+fup(3)*fup(2)*flo(1)&
                        *in%img_r8(ilo(1),iup(2)+1,iup(3)+1)
                   do ixx=ilo(1)+1,iup(1)
                      s=s+fup(3)*fup(2)*in%img_r8(ixx,iup(2)+1,iup(3)+1)
                   end do
                   if(fup(1)/=0) &
                        s=s+fup(3)*fup(2)*fup(1)&
                        *in%img_r8(iup(1)+1,iup(2)+1,iup(3)+1)
                end if
             endif
             out%img_r8(ix,iy,iz)=s*norm
!!$#ifndef 8
!!$ out%img_r8(ix,iy,iz)=s*norm
!!$#else
!!$ out%img_r8(ix,iy,iz)=nint(s*norm)
!!$#endif
          end do
       end do
    end do
  end subroutine
  !------------------------------------------------------!
  subroutine mv_img_r8(rtp,int_type,in,out)
    ! Transformation is defined at mv_img.
    type(rigid_transform_param),intent(in)::rtp
    integer*4,intent(in)::int_type
    type(image_data),intent(in)::in
    type(image_data),intent(out)::out
    integer*4 ix,iy,iz
    real*8 a(3,4)
    real*8 xin(3),xout(3)
    call alloc_image_data(out,out%n,(64))
    call set_affine(rtp,a)
!!$ call set_affine_rigid(rtp,a)
    do iz=1,out%n(3)
       xout(3)=iz*out%d(3)-out%offset(3)
       do iy=1,out%n(2)
          xout(2)=iy*out%d(2)-out%offset(2)
          do ix=1,out%n(1)
             xout(1)=ix*out%d(1)-out%offset(1)
             xin=((a.by.xout)+in%offset)/in%d
             select case(int_type)
             case(MV_INT_TYPE_NN)
                out%img_r8(ix,iy,iz)=&
                     & get_NN_interpolated_r8(in%n,xin,in%img_r8)
             case(MV_INT_TYPE_L2p)
                out%img_r8(ix,iy,iz)=&
                     & get_L2p_interpolated_r8(in%n,xin,in%img_r8)
             case(MV_INT_TYPE_SINC)
                out%img_r8(ix,iy,iz)=&
                     & get_sinc_interpolated_r8(in%n,xin,in%img_r8)
             case default
                out%img_r8(ix,iy,iz)=&
                     & get_Lnp_interpolated_r8(in%n,xin,in%img_r8)
             end select
          end do
       end do
    end do
  end subroutine
  !------------------------------------------------------!
  subroutine kmclusterize_img_r8(k,pmin,pmax,errrel_min,n,img,cimg)
    integer*4,intent(in)::k ! number of cluster
    real*8,intent(in)::pmin,pmax
    real*8,intent(in)::errrel_min
    integer*4,intent(in)::n(3)
    real(8),intent(in)::img(:,:,:)
    integer*4,intent(out)::cimg(:,:,:)
    real*8,allocatable::cmean(:)
    integer*4,parameter::NUM_ITERATION_MAX=100
    integer*4 nh
    real*8,allocatable::h(:)
    logical*4,allocatable::cstat(:)
    real*8 x,xx,dx
    integer*4 i
    nh=idnint(pmax-pmin)
    allocate(h(nh))
    call get_hist_r8(n,img,pmin,pmax,nh,h)
    allocate(cmean(k))
    allocate(cstat(k))
    dx=dble(nh)/k
    do i=1,k
       cmean(i)=i*dx-dx/2
    end do
    cstat=.true.
    xx=1
    i=0
    do i=1,NUM_ITERATION_MAX
       call set_cluster_mean
       call kmclustering_cost(x)
       if(dabs(x/xx-1)<errrel_min) exit
       xx=x
    end do
    call kmclassify()
    deallocate(h,cmean,cstat)
  contains
    integer*4 function get_upper_boundary(ipar0)
      integer*4,intent(in)::ipar0
      integer*4 ipar
      get_upper_boundary=nh
      do ipar=ipar0+1,k
         if(cstat(ipar)) then
            get_upper_boundary=max(1,min(nh,floor((cmean(ipar0)+cmean(ipar))/2)))
            exit
         end if
      end do
    end function get_upper_boundary
    !
    subroutine set_cluster_mean
      integer*4 ih,ipar
      integer*4 lower,upper,s1,s2
      lower=1
      do ipar=1,k
         if(.not.cstat(ipar)) cycle
         upper=get_upper_boundary(ipar)
         s1=0
         s2=0
         do ih=lower,upper
            s1=s1+h(ih)*ih
            s2=s2+h(ih)
         end do
         if(s2/=0) then
            cmean(ipar)=s1/s2
         else
            cstat(ipar)=.false.
         endif
         if(upper==nh) exit
         lower=upper+1
      end do
    end subroutine set_cluster_mean
    !
    subroutine kmclustering_cost(kmc)
      real*8,intent(out)::kmc
      integer*4 ipar,lower,upper,ih
      kmc=0
      lower=1
      do ipar=1,k
         if(.not.cstat(ipar)) cycle
         upper=get_upper_boundary(ipar)
         do ih=lower,upper
            kmc=kmc+h(ih)*(ih-cmean(ipar))**2
         end do
         if(upper==nh) exit
         lower=upper+1
      end do
    end subroutine kmclustering_cost
    !
    subroutine kmclassify()
      real*8 pval,d,dmin
      integer*4 icmin,ic
      integer*4 ip
      do ip=1,n(1)*n(2)*n(3)
         pval=img(ip,1,1)
         if(pval>=pmax.or.pval<=pmin)then
            cimg(ip,1,1)=0
            cycle
         endif
         dmin=HUGE(0.D0)
         icmin=0
         do ic=1,k
            if(.not.cstat(ic)) cycle
            d=dabs(cmean(ic)-pval)
            if(d<dmin) then
               dmin=d
               icmin=ic
            endif
            if(pval<cmean(ic)) exit
         end do
         cimg(ip,1,1)=icmin
      end do
    end subroutine kmclassify
  end subroutine
  !===========================================================================!
  subroutine contrast_norm_img_r8(nx,ny,nz,lo,up,n,img)
    ! contrast equalization
    integer*4,intent(in)::nx,ny,nz
    real*8,intent(in)::lo,up
    integer*4,intent(in)::n(3)
    real(8),intent(inout)::img(:,:,:)
    real*8,allocatable::mimg(:,:,:)
    integer*4 ix,iy,iz,k
    integer*4 iix,iiy,iiz
    integer*4 iix1,iix2,iiy1,iiy2,iiz1,iiz2
    real*8 s,s2,x
    allocate(mimg(n(1),n(2),n(3)))
    do iz=1,n(3)
       iiz1=max(iz-nz,1)
       iiz2=min(iz+nz,n(3))
       do iy=1,n(2)
          iiy1=max(iy-ny,1)
          iiy2=min(iy+ny,n(2))
          do ix=1,n(1)
             iix1=max(ix-nx,1)
             iix2=min(ix+nx,n(1))
             k=0
             s=0
             do iiz=iiz1,iiz2
                do iiy=iiy1,iiy2
                   do iix=iix1,iix2
                      x=img(iix,iiy,iiz)
                      if(x<lo.or.x>up) cycle
                      k=k+1
                      s=s+x
                   end do
                end do
             end do
             if(k/=0) then
                mimg(ix,iy,iz)=s/dble(k)
             else
                mimg(ix,iy,iz)=lo
             endif
          end do
       end do
    end do
    do iz=1,n(3)
       iiz1=max(iz-nz,1)
       iiz2=min(iz+nz,n(3))
       do iy=1,n(2)
          iiy1=max(iy-ny,1)
          iiy2=min(iy+ny,n(2))
          do ix=1,n(1)
             iix1=max(ix-nx,1)
             iix2=min(ix+nx,n(1))
             k=0
             s=0
             s2=0
             do iiz=iiz1,iiz2
                do iiy=iiy1,iiy2
                   do iix=iix1,iix2
                      x=mimg(iix,iiy,iiz)
                      if(x<lo.or.x>up) cycle
                      k=k+1
                      s=s+x
                      s2=s2+x**2
                   end do
                end do
             end do
             if(k/=0) then
                img(ix,iy,iz)=mimg(ix,iy,iz)&
                     /dsqrt(s2/dble(k)-(s/dble(k))**2)
             else
                img(ix,iy,iz)=lo
             endif
          end do
       end do
    end do
  end subroutine
  !
  real*8 function get_chi2_w_affine_r8(a,tar,ref) ! vulnerable
    ! to be a object function for mi-base matching problem
    real*8,intent(in)::a(3,4)
    type(image_data),intent(in)::tar
    type(image_data),intent(in)::ref
    integer*4,parameter::IREF=2,ITAR=1
    integer*4 ix,iy,iz
    real*8 xtar(3),xref(3)
    real*8 p(2)
    real*8 chisqr1,chisqr2
    integer*4,pointer::int_type
    int_type=>bws_gmwrd%ipd%it
    chisqr1=0
    chisqr2=0
    do iz=1,ref%n(3)
       xref(3)=iz*ref%d(3)-ref%offset(3)
       do iy=1,ref%n(2)
          xref(2)=iy*ref%d(2)-ref%offset(2)
          do ix=1,ref%n(1)
             xref(1)=ix*ref%d(1)-ref%offset(1)
             p(IREF)=ref%img_r8(ix,iy,iz)
             xtar=((a.by.xref)+tar%offset)/tar%d
             select case(int_type)
             case(MV_INT_TYPE_L2p)
                p(ITAR)=get_L2p_interpolated_r8(tar%n,xtar,tar%img_r8)
             case(MV_INT_TYPE_NN)
                p(ITAR)=get_NN_interpolated_r8(tar%n,xtar,tar%img_r8)
             case(MV_INT_TYPE_SINC)
                p(ITAR)=get_SINC_interpolated_r8(tar%n,xtar,tar%img_r8)
             case default
                p(ITAR)=get_Lnp_interpolated_r8(tar%n,xtar,tar%img_r8)
             end select
             chisqr1=chisqr1+(p(IREF)-p(ITAR))**2
          end do
       end do
    end do
    !
    do iz=1,tar%n(3)
       xtar(3)=iz*tar%d(3)-tar%offset(3)
       do iy=1,tar%n(2)
          xtar(2)=iy*tar%d(2)-tar%offset(2)
          do ix=1,tar%n(1)
             xtar(1)=ix*tar%d(1)-tar%offset(1)
             p(ITAR)=tar%img_r8(ix,iy,iz)
             xref=((a.invby.xtar)+ref%offset)/ref%d
             select case(int_type)
             case(MV_INT_TYPE_L2p)
                p(IREF)=get_L2p_interpolated_r8(ref%n,xref,ref%img_r8)
             case(MV_INT_TYPE_NN)
                p(IREF)=get_NN_interpolated_r8(ref%n,xref,ref%img_r8)
             case(MV_INT_TYPE_SINC)
                p(IREF)=get_SINC_interpolated_r8(ref%n,xref,ref%img_r8)
             case default
                p(IREF)=get_Lnp_interpolated_r8(ref%n,xref,ref%img_r8)
             end select
             chisqr2=chisqr2+(p(IREF)-p(ITAR))**2
          end do
       end do
    end do
    get_chi2_w_affine_r8=chisqr1+chisqr2
  end function get_chi2_w_affine_r8
  !
  real*8 function get_H_r8(n,img,ns,a,b)
    integer*4,intent(in)::n(3)
    real(8) img(:,:,:)
    integer*4,intent(in)::ns
    real*8,intent(in)::a,b
    real*8,allocatable::p(:)
    integer*4 i,k
    real*8 H,s
    allocate(p(ns))
    p=0
    do i=1,product(n)
       k=p2i(img(i,1,1))
       p(k)=p(k)+1
    end do
    s=sum(p)
    if(s/=0) then
       p=p/sum(p)
       H=1
       do k=1,ns
          if(p(k)==0) cycle
          H=H*p(k)**p(k)
       end do
       get_H_r8=log(H)
    else
       get_H_r8=0
    end if
    deallocate(p)
  contains
    integer*4 function p2i(x)
      real(8),intent(in)::x
      p2i=idnint(x*a+b)
      if(p2i<=0) then
         p2i=1
      elseif(p2i>ns)then
         p2i=ns
      endif
    end function p2i
    !
  end function get_H_r8
  !
  real*8 function get_mi_w_affine_r8(a,tar,ref) ! vulnerable
    ! to be a object function for mi-base matching problem
    real*8,intent(in)::a(3,4)
    type(image_data),intent(in)::tar
    type(image_data),intent(in)::ref
    integer*4,parameter::IREF=2,ITAR=1
    integer*4 ix,iy,iz
    real*8 xtar(3),xref(3)
    real*8 p(2)
    real*8,pointer::h(:,:),hh(:,:)
    integer*4,pointer::int_type,flag
    integer*4 k,nns(3,32)
    real*8 ws(32)
    integer*4 iip(2)
    real*8 sh
    h=>bws_h(:,:,1) ! defined at butyD
    hh=>bws_h(:,:,2)
    int_type=>bws_gmwrd%ipd%it
    flag=>bws_gmwrd%flag
    h=0
    do iz=1,ref%n(3)
       xref(3)=iz*ref%d(3)-ref%offset(3)
       do iy=1,ref%n(2)
          xref(2)=iy*ref%d(2)-ref%offset(2)
          do ix=1,ref%n(1)
             xref(1)=ix*ref%d(1)-ref%offset(1)
             xtar=((a.by.mkrand(xref,ref%d(1)))+tar%offset)/tar%d
             xtar=mkrand(xtar,tar%d(1))
             if(.not.IsInside(xtar,tar%n)) then

                iip(IREF)=adj_p(ref%img_r8(ix,iy,iz),ref%a,ref%b,ref%ns)



                h(1,iip(IREF))=h(1,iip(IREF))+1
             else
                if(int_type==HIST_INT_TYPE_PV) then

                   iip(IREF)=adj_p(ref%img_r8(ix,iy,iz),ref%a,ref%b,ref%ns)



                   do k=1,get_L2P_weight(xtar,tar%n)
                      if(ws(k)==0) cycle

                      iip(ITAR)=adj_p(tar%img_r8(nns(1,k),nns(2,k),nns(3,k)),&
                           tar%a,tar%b,tar%ns)




                      h(iip(ITAR),iip(IREF))=h(iip(ITAR),iip(IREF))+ws(k)
                   end do
                else
                   p(IREF)=ref%img_r8(ix,iy,iz)
                   p(ITAR)=get_intp(xtar,tar)
                   call histograming
                end if
             end if
          end do
       end do
    end do
    h=>bws_h(:,:,2)
    h=0
    do iz=1,tar%n(3)
       xtar(3)=iz*tar%d(3)-tar%offset(3)
       do iy=1,tar%n(2)
          xtar(2)=iy*tar%d(2)-tar%offset(2)
          do ix=1,tar%n(1)
             xtar(1)=ix*tar%d(1)-tar%offset(1)
             xref=((a.invby.mkrand(xtar,tar%d(1)))+ref%offset)/ref%d
             xref=mkrand(xref,ref%d(1))
             if(.not.IsInside(xref,ref%n)) then

                iip(ITAR)=adj_p(tar%img_r8(ix,iy,iz),tar%a,tar%b,tar%ns)



                h(iip(ITAR),1)=h(iip(ITAR),1)+1
             else
                if(int_type==HIST_INT_TYPE_PV) then

                   iip(ITAR)=adj_p(tar%img_r8(ix,iy,iz),tar%a,tar%b,tar%ns)



                   do k=1,get_L2P_weight(xref,ref%n)
                      if(ws(k)==0) cycle

                      iip(IREF)=adj_p(ref%img_r8(nns(1,k),nns(2,k),nns(3,k)),&
                           ref%a,ref%b,ref%ns)




                      h(iip(ITAR),iip(IREF))=h(iip(ITAR),iip(IREF))+ws(k)
                   end do
                else
                   p(ITAR)=tar%img_r8(ix,iy,iz)
                   p(IREF)=get_intp(xref,ref)
                   call histograming
                end if
             end if
          end do
       end do
    end do
    h=>bws_h(:,:,1)
    h=h+hh
    sh=sum(h)
    if(sh/=0) then
       h=h/sh
       if(iand(flag,GMWRF_USE_SMOOTH_PDF)/=0) then
          call smooth2d(tar%ns,ref%ns,h,hh)
          h=>hh
       endif
       get_mi_w_affine_r8=eval_mutual_information(h)
    else
       get_mi_w_affine_r8=0
    end if
    !
  contains
    !
    function mkrand(x0,r)
      real*8 mkrand(3)
      real*8,intent(inout)::x0(3)
      real*8,intent(in)::r
      real*8 ran(4)
      mkrand=x0
      return
!!$ call random_number(ran)
!!$ ran(1)=sqrt(-2*log(1-ran(1)))*cos(pid*ran(4))
!!$ mkrand=x0+r*(ran-.5)/4.
!!$ return
      ran(1)=ran(1)*r/8.
      ran(2)=(ran(2)-.5)*2.*pi
      ran(3)=ran(3)*pi
      mkrand(1)=x0(1)+ran(1)*sin(ran(3))*cos(ran(2))
      mkrand(2)=x0(2)+ran(1)*sin(ran(3))*sin(ran(2))
      mkrand(3)=x0(3)+ran(1)*cos(ran(3))
    end function mkrand
    !
    real*8 function get_intp(x0,id0)
      real*8,intent(in)::x0(3)
      type(image_data),intent(in)::id0
      select case(int_type)
      case(MV_INT_TYPE_L2p)
         get_intp=get_L2p_interpolated_r8(id0%n,x0,id0%img_r8)
      case(MV_INT_TYPE_NN)
         get_intp=get_NN_interpolated_r8(id0%n,x0,id0%img_r8)
      case(MV_INT_TYPE_SINC)
         get_intp=get_SINC_interpolated_r8(id0%n,x0,id0%img_r8)
      case default
         get_intp=get_Lnp_interpolated_r8(id0%n,x0,id0%img_r8)
      end select
    end function get_intp
    !
    logical*4 function isInside(x0,n0)
      real*8,intent(in)::x0(3)
      integer*4,intent(in)::n0(3)
      integer*4 j
      isInside=.false.
      do j=1,3
         if(x0(j)<=0) then
            return
         elseif(idnint(x0(j))<=0.or.idnint(x0(j))>n0(j)) then
            return
         endif
      end do
      isInside=.true.
    end function isInside
    !
    integer*4 function get_L2P_weight(x0,n)
      real*8,intent(in)::x0(3)
      integer*4,intent(in)::n(3)
      integer*4 ix0(3)
      integer*4 skip,istat(3)
      integer*4 j,k
      real*8 fx(3),qx(3)
      get_L2P_weight=0
      ix0=int(x0)
      skip=0
      ws=0
      k=0
      istat=0
      nns=0
      !
      do j=1,3
         if(x0(j)<0)then
            skip=1
            return
         elseif(ix0(j)>0.and.ix0(j)<n(j)) then
            istat(j)=0
         elseif(ix0(j)<0.or.ix0(j)>n(j)) then
            skip=1
            return
         elseif(ix0(j)==n(j)) then
            istat(j)=2
         else
            istat(j)=1
         endif
      end do
      fx=x0-dble(ix0)
      qx=1-fx
      if(iand(ior(istat(3),ior(istat(1),istat(2))),1)==0) then
         k=k+1
         nns(:,k)=ix0
         ws(k)=qx(1)*qx(2)*qx(3)
      end if
      if(iand(istat(1),2)==0.and.iand(ior(istat(2),istat(3)),1)==0&
           .and.fx(1)/=0) then
         k=k+1
         nns(1,k)=ix0(1)+1
         nns(2,k)=ix0(2)
         nns(3,k)=ix0(3)
         ws(k)=fx(1)*qx(2)*qx(3)
      end if
      if(iand(istat(2),2)==0.and.iand(ior(istat(1),istat(3)),1)==0&
           .and.fx(2)/=0) then
         k=k+1
         nns(1,k)=ix0(1)
         nns(2,k)=ix0(2)+1
         nns(3,k)=ix0(3)
         ws(k)=qx(1)*fx(2)*qx(3)
      end if
      if(iand(istat(3),1)==0.and.iand(ior(istat(1),istat(2)),2)==0&
           .and.fx(2)/=0.and.fx(1)/=0) then
         k=k+1
         nns(1,k)=ix0(1)+1
         nns(2,k)=ix0(2)+1
         nns(3,k)=ix0(3)
         ws(k)=fx(1)*fx(2)*qx(3)
      end if
      !
      if(iand(istat(3),2)==0.and.iand(ior(istat(1),istat(2)),1)==0&
           .and.fx(3)/=0) then
         k=k+1
         nns(1,k)=ix0(1)
         nns(2,k)=ix0(2)
         nns(3,k)=ix0(3)+1
         ws(k)=qx(1)*qx(2)*fx(3)
      end if
      !
      if(iand(istat(1),1)==0.and.iand(ior(istat(2),istat(3)),2)==0&
           .and.fx(2)/=0.and.fx(3)/=0) then
         k=k+1
         nns(1,k)=ix0(1)
         nns(2,k)=ix0(2)+1
         nns(3,k)=ix0(3)+1
         ws(k)=qx(1)*fx(2)*fx(3)
      end if
      !
      if(iand(istat(2),1)==0.and.iand(ior(istat(1),istat(3)),2)==0&
           .and.fx(1)/=0.and.fx(3)/=0) then
         k=k+1
         nns(1,k)=ix0(1)+1
         nns(2,k)=ix0(2)
         nns(3,k)=ix0(3)+1
         ws(k)=fx(1)*qx(2)*fx(3)
      end if
      !
      if(iand(ior(istat(3),ior(istat(1),istat(2))),2)==0&
           .and.fx(1)/=0.and.fx(2)/=0.and.fx(3)/=0) then
         k=k+1
         nns(:,k)=ix0+1
         ws(k)=fx(1)*fx(2)*fx(3)
      end if
      !
!!$ s=sum(ws)
!!$ if(s/=0) then
!!$ ws=ws/s
      get_L2P_weight=k
!!$ end if
    end function get_L2P_weight
    !
    !
    integer*4 function adj_p(p0,a0,b0,n0)
      real*8,intent(in)::p0
      real*8,intent(in)::a0,b0
      integer*4,intent(in)::n0
      adj_p=idnint(p0*a0+b0)
      if(adj_p<=0) then
         adj_p=1
      elseif(adj_p>n0)then
         adj_p=n0
      endif
    end function adj_p
    !
    subroutine histograming
      real*8 fp(2),qp(2)
      integer*4 ip(2)
      integer*4 hstat
      ip(ITAR)=adj_p(p(ITAR),tar%a,tar%b,tar%ns)
      ip(IREF)=adj_p(p(IREF),ref%a,ref%b,ref%ns)
      h(ip(ITAR),ip(IREF))=h(ip(ITAR),ip(IREF))+1
      ip=int(p)
      fp=p-ip
      hstat=0
      if(ip(ITAR)<=0) then
         ip(ITAR)=1
         fp(ITAR)=0
      end if
      if(ip(ITAR)>=tar%ns) then
         ip(ITAR)=tar%ns
         fp(ITAR)=0
         hstat=1
      end if
      !
      if(ip(IREF)<=0) then
         ip(IREF)=1
         fp(IREF)=0
      end if
      if(ip(IREF)>=ref%ns) then
         ip(IREF)=ref%ns
         fp(IREF)=0
         hstat=ior(hstat,2)
      end if
      !
      qp=1-fp
      h(ip(ITAR),ip(IREF))=h(ip(ITAR),ip(IREF))+qp(ITAR)*qp(IREF)
      select case(hstat)
      case(0)
         h(ip(ITAR)+1,ip(IREF))=h(ip(ITAR)+1,ip(IREF))+fp(ITAR)*qp(IREF)
         h(ip(ITAR),ip(IREF)+1)=h(ip(ITAR),ip(IREF)+1)+qp(ITAR)*fp(IREF)
         h(ip(ITAR)+1,ip(IREF)+1)=h(ip(ITAR),ip(IREF)+1)+fp(ITAR)*fp(IREF)
      case(1)
         h(ip(ITAR),ip(IREF)+1)=h(ip(ITAR),ip(IREF)+1)+qp(ITAR)*fp(IREF)
      case(2)
         h(ip(ITAR)+1,ip(IREF))=h(ip(ITAR)+1,ip(IREF))+fp(ITAR)*qp(IREF)
      end select
    end subroutine histograming
    !-=====================================================
    real*8 function eval_mutual_information(pdf)
      real*8,intent(in)::pdf(:,:)
      real*8 H1,H2,H12
      real*8 p1,p2,p12
      integer*4 it,ir
      H1=1
      do ir=1,ref%ns
         p1=sum(pdf(:,ir))
         if(p1==0) cycle
         H1=H1*p1**p1
      end do
      H2=1
      do it=1,tar%ns
         p2=sum(pdf(it,:))
         if(p2==0) cycle
         H2=H2*p2**p2
      end do
      H12=1
      do ir=1,ref%ns
         do it=1,tar%ns
            p12=pdf(it,ir)
            if(p12==0) cycle
            H12=H12*p12**p12
         end do
      end do
      if(H12/=0) then
         eval_mutual_information=-log(H1*H2)/log(H12)
      else
         eval_mutual_information=0
      endif
    end function eval_mutual_information
    !
  end function get_mi_w_affine_r8
  !
  subroutine infopot_img_r8(s,ns,n,img,buf)
    ! ns must be >= up-lo
    use muty
    ! evaluate information potential of an image
    real*8,intent(in)::s
    integer*4,intent(in)::ns
    integer*4,intent(in)::n(3)
    real(8),intent(in)::img(:,:,:)
    real(8),intent(out)::buf(:,:,:)
    integer*4 ix,iy,iz
    integer*4 i,j
    real*8 sigmasqr,val,norm,diag
    integer*4,parameter::num_sample=8
    integer*4,parameter::n_tbl=8
    real*8,parameter::dt_tbl=1/dble(n_tbl)
    real(8) xij(num_sample)
    real*8,allocatable::gi(:)
    sigmasqr=2*s**2
    norm=1/dsqrt(2*PI*sigmasqr)
    diag=num_sample*norm
    norm=norm*2
    allocate(gi(0:ns))
    do i=0,ns
       gi(i)=g(dble(i)*dt_tbl)*norm ! expects integer image
    end do
    do iz=1,n(3)-1
       do iy=1,n(2)-1
          do ix=1,n(1)-1
             xij(1:2)=img(ix:ix+1,iy,iz)
             xij(3:4)=img(ix:ix+1,iy+1,iz)
             xij(5:6)=img(ix:ix+1,iy,iz+1)
             xij(7:8)=img(ix:ix+1,iy+1,iz+1)
             val=diag
             do i=1,8
                do j=i+1,8
                   val=val+gi(abs(n_tbl*int((xij(i)-xij(j)))))
                end do
             end do
             buf(ix,iy,iz)=val
          end do
       end do
    end do
    deallocate(gi)
  contains
    real*8 function g(xarg)
      real*8,intent(in)::xarg
      g=dexp(-xarg**2/(2*sigmasqr))
    end function g
  end subroutine
  !
  real(8) function get_NN_interpolated_r8(n,x,img)
    integer*4,intent(in)::n(3)
    real*8,intent(in)::x(3)
    real(8),intent(in)::img(:,:,:)
    integer*4 ix(3)
    integer*4 skip
    integer*4 j
    ix=idnint(x)
    skip=0
    do j=1,3
       if(x(j)<0) then
          skip=1
          exit
       elseif(ix(j)<=0.or.ix(j)>n(j))then
          skip=1
          exit
       endif
    end do
    if(skip==0) then
       get_NN_interpolated_r8=img(ix(1),ix(2),ix(3))
    else
       get_NN_interpolated_r8=PVAL_OUT_OF_SIGHT
    endif
  end function
  !
  real(8) function get_Lnp_interpolated_r8(n,x,img)
    ! 3D analogue of
    ! Milton Abramowitz and Irene A. Stegun (1972).
    ! "Lagrange interpolation, Equally Spaced Abscissas, n Point Formula"
    ! HMF. 25.2.6-7
    integer*4,intent(in)::n(3)
    real*8,intent(in)::x(3)
    real(8),intent(in)::img(:,:,:)
    integer*4,pointer:: mh(:),lo(:),up(:),tlo(:),tup(:)
    real*8,pointer::nu(:,:)
    real*8 fx(3)
    integer*4 ix(3)
    integer*4 jx,jy,jz,kx,ky,kz
    real*8 s,Az,Ay
    mh=>bws_lnpt%nh
    lo=>bws_lnpt%lo
    up=>bws_lnpt%nh
    tlo=>bws_lnpt%tlo
    tup=>bws_lnpt%tup
    nu=>bws_lnpt%nu
    s=0
    ix=int(x)
    fx=x-dble(ix)
    do kz=lo(3),up(3)
       jz=ix(3)+kz
       if(jz<=0) cycle
       if(jz>n(3)) exit
       Az=Akn(kz,3,fx(3))
       if(Az==0) cycle
       do ky=lo(2),up(2)
          jy=ix(2)+ky
          if(jy<=0) cycle
          if(jy>n(2)) exit
          Ay=Akn(ky,2,fx(2))*Az
          if(Ay==0) cycle
          do kx=lo(1),up(1)
             jx=ix(1)+kx
             if(jx<=0) cycle
             if(jx>n(1)) exit
             s=s+img(jx,jy,jz)*Akn(kx,1,fx(1))*Ay
          end do
       end do
    end do
    get_Lnp_interpolated_r8=s
  contains
    real*8 function Akn(k,in,p)
      integer*4,intent(in)::k
      integer*4,intent(in)::in
      real*8,intent(in)::p
      integer*4 t
      Akn=nu(k,in)
      do t=tlo(in),tup(in)
         if(mh(in)-t==-k) cycle
         Akn=Akn*(p+mh(in)-t)
      end do
    end function Akn
  end function
  !
  real(8) function get_L2P_interpolated_r8(n,x,img)
    ! 3D analogue of
    ! Milton Abramowitz and Irene A. Stegun (1972).
    ! "Bivariate Interpolation, Four Point Formula."
    ! HMF. 25.2.66
    integer*4,intent(in)::n(3)
    real*8,intent(in)::x(3)
    real(8),intent(in)::img(:,:,:)
    integer*4 j
    real*8 fx(3),qx(3),s
    integer*4 ix(3)
    integer*4 skip,istat(3)
    ix=int(x)
    skip=0
    do j=1,3
       if(x(j)<0)then
          skip=1
          exit
       elseif(ix(j)>0.and.ix(j)<n(j)) then
          istat(j)=0
       elseif(ix(j)<0.or.ix(j)>n(j)) then
          skip=1
          exit
       elseif(ix(j)==n(j)) then
          istat(j)=2
       else
          istat(j)=1
       endif
    end do
    if(skip==0) then
       fx=x-dble(ix)
       qx=1-fx
       s=0
       if(iand(ior(istat(3),ior(istat(1),istat(2))),1)==0) &
            s=s+img(ix(1),ix(2),ix(3)) &
            *qx(1)*qx(2)*qx(3)
       if(iand(istat(1),2)==0.and.iand(ior(istat(2),istat(3)),1)==0) &
            s=s+img(ix(1)+1,ix(2),ix(3)) &
            *fx(1)*qx(2)*qx(3)
       if(iand(istat(2),2)==0.and.iand(ior(istat(1),istat(3)),1)==0) &
            s=s+img(ix(1),ix(2)+1,ix(3)) &
            *qx(1)*fx(2)*qx(3)
       if(iand(istat(3),1)==0.and.iand(ior(istat(1),istat(2)),2)==0) &
            s=s+img(ix(1)+1,ix(2)+1,ix(3)) &
            *fx(1)*fx(2)*qx(3)
       if(iand(istat(3),2)==0.and.iand(ior(istat(1),istat(2)),1)==0) &
            s=s+img(ix(1),ix(2),ix(3)+1) &
            *qx(1)*qx(2)*fx(3)
       if(iand(istat(1),1)==0.and.iand(ior(istat(2),istat(3)),2)==0) &
            s=s+img(ix(1),ix(2)+1,ix(3)+1) &
            *qx(1)*fx(2)*fx(3)
       if(iand(istat(2),1)==0.and.iand(ior(istat(1),istat(3)),2)==0) &
            s=s+img(ix(1)+1,ix(2),ix(3)+1) &
            *fx(1)*qx(2)*fx(3)
       if(iand(ior(istat(3),ior(istat(1),istat(2))),2)==0) &
            s=s+img(ix(1)+1,ix(2)+1,ix(3)+1) &
            *fx(1)*fx(2)*fx(3)
       get_L2P_interpolated_r8=s
    else
       get_L2P_interpolated_r8=PVAL_OUT_OF_SIGHT
    end if
  end function
  !
  real(8) function get_sinc_interpolated_r8(n,x,img)
    integer*4,intent(in)::n(3)
    real*8,intent(in)::x(3)
    real(8),intent(in)::img(:,:,:)
    integer*4 jx,jy,jz
    integer*4 ix(3),iix(3),lo(3),up(3)
    real*8 fx(3)
    real*8 s
    real*8 sk(3)
    integer*4,pointer::d_sinc_kernel_tbl,n_sinc_mask(:)
    real*8,pointer::sinc_kernel_tbl(:)
    d_sinc_kernel_tbl => bws_sinct%d
    n_sinc_mask => bws_sinct%m
    sinc_kernel_tbl => bws_sinct%kt
    ix=int(x)
    fx=abs(x-ix)*d_sinc_kernel_tbl
    s=0
    lo=max(1,ix-n_sinc_mask)-ix
    up=min(n,ix+n_sinc_mask)-ix
    if(up(1)/=lo(1).and.up(2)/=lo(2).and.up(3)/=lo(3)) then
       do jz=lo(3),up(3)
          iix(3)=ix(3)+jz
          sk(3)=sinc_kernel_tbl(abs(idnint(fx(3))-jz*d_sinc_kernel_tbl))
          do jy=lo(2),up(2)
             iix(2)=ix(2)+jy
             sk(2)=sinc_kernel_tbl(abs(idnint(fx(2))-jy*d_sinc_kernel_tbl))&
                  *sk(3)
             do jx=lo(1),up(1)
                iix(1)=ix(1)+jx
                sk(1)=sinc_kernel_tbl(abs(idnint(fx(1))-jx*d_sinc_kernel_tbl))
                s=s+img(iix(1),iix(2),iix(3))*sk(1)*sk(2)
             end do
          end do
       end do
       get_sinc_interpolated_r8=s
    else
       get_sinc_interpolated_r8=PVAL_OUT_OF_SIGHT
    endif
  end function
  !
  subroutine adapbin_filter_r8(n,img,mask)
    ! adaptive binary-zation
    integer*4,intent(in)::n(3)
    real(8),intent(inout)::img(:,:,:)
    real(8),allocatable::buf(:,:,:)
    integer*4,intent(in)::mask(3)
    integer*4 ix,iy,iz,jx,jy,jz,jx1,jx2,jy1,jy2,jz1,jz2
    integer*4 my,mz
    real*8 norm
    allocate(buf(n(1),n(2),n(3)))
    buf=0
    do iz=1,n(3)
       jz1=max(1,iz-mask(3))
       jz2=min(n(3),iz+mask(3))
       mz=jz2-jz1+1
       do iy=1,n(2)
          jy1=max(1,iy-mask(2))
          jy2=min(n(2),iy+mask(2))
          my=mz*(jy2-jy1+1)
          do ix=1,n(1)
             jx1=max(1,ix-mask(1))
             jx2=min(n(1),ix+mask(1))
             norm=1/dble((jx2-jx1+1)*my)
             do jz=jz1,jz2
                do jy=jy1,jy2
                   do jx=jx1,jx2
                      buf(ix,iy,iz)=buf(ix,iy,iz)+img(jx,jy,jz)
                   end do
                end do
             end do
             buf(ix,iy,iz)=buf(ix,iy,iz)*norm
          end do
       end do
    end do
    do iz=1,n(3)
       do iy=1,n(2)
          do ix=1,n(1)
             if(buf(ix,iy,iz)>img(ix,iy,iz)) then
                img(ix,iy,iz)=0
             else
                img(ix,iy,iz)=1
             end if
          end do
       end do
    end do
    deallocate(buf)
  end subroutine
 !
  subroutine get_cm_img_r8(n,d,o,img,cm)
    integer*4,intent(in)::n(3)
    real*8,intent(in)::d(3)
    real*8,intent(in)::o(3)
    real(8),intent(in)::img(:,:,:)
    real*8,intent(out)::cm(3)
    integer*4 ix,iy,iz
    real*8 x(3),m
    m=0
    cm=0
    do iz=1,n(3)
       x(3)=iz*d(3)-o(3)
       do iy=1,n(2)
          x(2)=iy*d(2)-o(2)
          do ix=1,n(1)
             x(1)=ix*d(1)-o(1)
             m=m+img(ix,iy,iz)
             cm=cm+img(ix,iy,iz)*x
          end do
       end do
    end do
    if(m/=0)then
       cm=cm/m
    else
       cm=0
    endif
  end subroutine
  !
  subroutine get_I_img_r8(n,d,o,cm,img,I)
    ! get moment of inertia
    integer*4,intent(in)::n(3)
    real*8,intent(in)::d(3)
    real*8,intent(in)::o(3)
    real*8,intent(in)::cm(3)
    real(8),intent(in)::img(:,:,:)
    real*8,intent(out)::I(3)
    integer*4 ix,iy,iz
    real*8 x(3),m
    m=0
    I=0
    do iz=1,n(3)
       x(3)=iz*d(3)-o(3)
       do iy=1,n(2)
          x(2)=iy*d(2)-o(2)
          do ix=1,n(1)
             x(1)=ix*d(1)-o(1)
             m=m+img(ix,iy,iz)
             I=I+img(ix,iy,iz)*(cm-x)**2
          end do
       end do
    end do
    if(m/=0)then
       I=I/m
    else
       I=0
    endif
  end subroutine
  !
  subroutine label2d_img_r8(in,out,maxkmax)
    type(image_data),intent(in)::in
    type(image_data),intent(out),target::out
    integer*4,intent(out)::maxkmax
    integer*4,parameter::MAX_LABEL=1024
    integer*4 ix,iy,iz
    integer*4 k,L,m,kmax
    integer*4,allocatable,target::label(:,:)
    integer*4,allocatable::eq(:)
    integer*4,pointer::Lp
    integer*4,pointer::out_img(:,:)
    call init_image_data(out)
    call cp_image_data(in,out)
    call alloc_image_data(out,in%n,DT_IKIND_4)
    allocate(label(in%n(1),in%n(2)))
    allocate(eq(0:MAX_LABEL))
    maxkmax=0
    do iz=1,in%n(3)
       eq=0
       label=0
       k=0
       nullify(Lp)
       out_img=>out%img_i4(1:in%n(1),1:in%n(2),iz)
       do iy=1,in%n(2)
          do ix=1,in%n(1)
             if(isValid(in%img_r8(ix,iy,iz))) then
                m=0
                L=0
                if(iy/=1) then
                   if(isValid(in%img_r8(ix,iy-1,iz))) then
                      m=m+1
                      Lp=>label(ix,iy-1)
                      L=Lp
                   end if
                   if(ix/=1) then
                      if(isValid(in%img_r8(ix-1,iy,iz))) then
                         Lp=>label(ix-1,iy)
                         if(m/=0.and.L/=Lp) call label_eq(L,Lp)
                         m=m+1
                         L=Lp
                      end if
                      if(isValid(in%img_r8(ix-1,iy-1,iz))) then
                         Lp=>label(ix-1,iy-1)
                         if(m/=0.and.L/=Lp) call label_eq(L,Lp)
                         m=m+1
                         L=Lp
                      end if
                   end if
                   if(ix/=in%n(1)) then
                      if(isValid(in%img_r8(ix+1,iy-1,iz))) then
                         Lp=>label(ix+1,iy-1)
                         if(m/=0.and.L/=Lp) call label_eq(L,Lp)
                         m=m+1
                         L=Lp
                      end if
                   end if
                else
                   if(ix/=1) then
                      if(isValid(in%img_r8(ix-1,iy,iz))) then
                         Lp=>label(ix-1,iy)
                         if(m/=0.and.L/=Lp) call label_eq(L,Lp)
                         m=m+1
                         L=Lp
                      end if
                   end if
                end if
                Lp=>label(ix,iy)
                if(m==0) then
                   if(k==MAX_LABEL) then
                      stop "MAX_LABEL"
                   endif
                   k=k+1
                   Lp=k
                else
                   Lp=L
                end if
             end if
          end do
       end do
       ! sorting equivalent labels
       kmax=kill_eq(k)
       maxkmax=max(kmax,maxkmax)
       ! update labels
       do iy=1,in%n(2)
          do ix=1,in%n(1)
             Lp=>label(ix,iy)
             Lp=eq(Lp)
          end do
       end do
       out_img=label
    end do
    deallocate(eq,label)
  contains
    logical*4 function isValid(p)
      real(8),intent(in)::p
      if(p>in%plo.and.p<in%pup) then
         isValid=.true.
      else
         isValid=.false.
      endif
    end function isValid
    !
    subroutine label_eq(L1,L2)
      integer*4,intent(in),target::L1,L2
      integer*4,pointer::el1,el2
      integer*4 ii,jj
      if(L1<=L2) then
         el1=>L1
         el2=>L2
      else
         el1=>L2
         el2=>L1
      end if
      jj=el2
      do
         if(eq(jj)==0) then
            if(eq(el1)==jj) return
            exit
         endif
         if(el1==eq(jj)) exit
         if(jj==eq(jj)) return
         if(jj<eq(jj)) then
            jj=eq(jj)
         else
            el2=eq(jj)
            if(el2<el1) then
               ii=el2
               el2=el1
               el1=ii
            end if
            jj=el2
         end if
      end do
      eq(jj)=el1
    end subroutine label_eq
    !
    integer*4 function kill_eq(kkmax)
      integer*4,intent(in)::kkmax
      integer*4 ii,jj,kk
      kk=0
      do ii=1,kkmax
         if(eq(ii)/=0) then
            jj=eq(ii)
            eq(ii)=eq(jj)
         else
            kk=kk+1
            eq(ii)=kk
         end if
      end do
      kill_eq=kk
    end function kill_eq
  end subroutine
 !
# 1711 "butyM.f90"
  !
  subroutine filter_img_r8(n,m,kern,img)
    integer*4,intent(in)::n(3) ! dimension(img)
    integer*4,intent(in)::m(3) ! kern(-m(1):m(1),...)
    real*8,intent(in)::kern(-m(1):m(1),-m(2):m(2),-m(3):m(3))
    real(8),intent(inout)::img(:,:,:)
    real(8),allocatable::out(:,:,:)
    integer*4 ix,iy,iz
    integer*4 ilo(3),iup(3)
    integer*4 jx,jy,jz,jjy,jjz,jjx
    real*8 s,cmax
    allocate(out(n(1),n(2),n(3)))
    cmax=(m(1)*2+1)*(m(2)*2+1)*(m(3)*2+1)
    ilo=1+m
    iup=n-m
    do iz=ilo(3),iup(3)
       do iy=ilo(2),iup(2)
          do ix=ilo(1),iup(1)
             s=0
             do jz=-m(3),m(3)
                jjz=iz+jz
                do jy=-m(2),m(2)
                   jjy=iy+jy
                   do jx=-m(1),m(1)
                      s=s+img(ix+jx,jjy,jjz)*kern(jx,jy,jz)
                   end do
                end do
             end do



             out(ix,iy,iz)=s

          end do
       end do
    end do
    ! bottom
    call conv(1,n(1), 1,n(2), 1,m(3))
    ! top
    call conv(1,n(1), 1,n(2), n(3)-m(3)+1,n(3))
    ! side y-1
    call conv(1,n(1), 1,m(2), m(3)+1,n(3)-m(3))
    ! side y+1
    call conv(1,n(1), n(2)-m(2)+1,n(2), m(3)+1,n(3)-m(3))
    ! side x-1
    call conv(1,m(1), m(2)+1,n(2)-m(2), m(3)+1,n(3)-m(3))
    ! side x+1
    call conv(n(1)-m(1)+1,n(1), m(2)+1,n(2)-m(2), m(3)+1,n(3)-m(3))
    img=out
    deallocate(out)
  contains
    subroutine conv(xlo,xup,ylo,yup,zlo,zup)
      integer*4,intent(in)::xlo,xup,ylo,yup,zlo,zup
      integer*4 jjx1,jjx2,jjy1,jjy2,jjz1,jjz2
      integer*4 c
      do iz=zlo,zup
         jjz1=max(1,iz-m(3))
         jjz2=min(n(3),iz+m(3))
         do iy=ylo,yup
            jjy1=max(1,iy-m(2))
            jjy2=min(n(2),iy+m(2))
            do ix=xlo,xup
               jjx1=max(1,ix-m(1))
               jjx2=min(n(1),ix+m(1))
               s=0
               c=0
               do jjz=jjz1,jjz2
                  do jjy=jjy1,jjy2
                     do jjx=jjx1,jjx2
                        s=s+img(jjx,jjy,jjz)*kern(jjx-ix,jjy-iy,jjz-iz)
                        c=c+1
                     end do
                  end do
               end do



               if(c/=0) out(ix,iy,iz)=s*cmax/c

            end do
         end do
      end do
    end subroutine conv
  end subroutine
  !
  subroutine rlcomp_img_r8(n,img)
    ! run-length compression of 0
    ! img must be positive or 0
    integer*4,intent(in)::n(3)
    real(8),intent(inout)::img(:,:,:)
    integer*4 i,j,m
    j=1
    m=0
    do i=1,product(n)
       if(img(i,1,1)<=0)then
          m=m+1
       elseif(m/=0) then
          img(j,1,1)=-m
          j=j+1
          m=0
       else
          img(j,1,1)=img(i,1,1)
          j=j+1
       endif
    end do
  end subroutine
  !
  subroutine rluncomp_img_r8(n,img,out)
    integer*4,intent(in)::n(3)
    real(8),intent(in)::img(:,:,:)
    real(8),intent(out)::out(:,:,:)
    integer*4 i,j,k
    j=0
    do i=1,product(n)
       if(img(i,1,1)<=0) then
          do k=1,int(abs(img(i,1,1)))
             j=j+1
             out(j,1,1)=0
          end do
       else
          j=j+1
          out(j,1,1)=img(i,1,1)
       end if
    end do
  end subroutine
end module
